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DEVELOPMENT  OF  PERTURBATION  PROCEDURES  FOR 
NONLINEAR  INVISCID  AND  VISCOUS  FLOWS 

by  David  Nixon 


ABSTRACT 

The  perturbation  theory  for  transonic  flow  is  further  devel¬ 
oped  for  two  and  three  dimensional  problems.  The  physical  per¬ 
turbation  theory,  in  which  the  perturbation  parameter  is  a 
physical  quantity,  is  applied  to  solutions  of  the  Navier-Stokes 
equations  in  two  dimensions.  The  mathematical  perturbation 
theory,  in  which  the  perturbation  parameter  is  a  measure  of  the 
difference  between  approximate  and  exact  solutions,  is  applied 
both  to  the  three-dimensional  potential  flow  problems  and  to  the 
two-dimensional  Navier-Stokes  equations.  The  strained  coordinate 
technique  is  used  to  treat  changes  in  location  of  any  shock  waves 
or  large  gradients.^ 

1 .  INTRODUCTION 


A  series  of  papers  recently  published  (refs.  1,  2,  and  3) 
concern  the  development  of  a  perturbation  method  for  transonic 
flow  problems.  A  basic  problem  in  transonic  flow  perturbations 
is  the  treatment  of  possible  movement  of  discontinuities  (shock 
waves)  which  can  make  the  perturbation  invalid.  A  means  of 
solving  this  problem,  using  a  strained  coordinate  method,  has 
been  derived  (ref.  1) .  The  location  of  the  discontinuities  in 
this  strained  coordinate  system  remain  the  same  throughout  the 
perturbation,  leading  to  a  final  linear  perturbation  equation 
for  the  perturbed  unknowns.  The  final  solution  is  not  linear 
because  of  the  transformation  from  the  strained  coordinate  system 
to  the  physical  coordinate  system.  However,  since  the  perturba¬ 
tion  equation  is  linear,  the  principle  of  superposition  can  be 
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used,  and  the  effect  of  several  different  types  of  perturbations 
can  be  easily,  and  inexpensively,  estimated. 

The  most  common  application  (ref.  2)  of  the  method  has  been 
concerned  with  extrapolating  transonic  flow  solutions.  For 
example,  given  two  transonic  flow  solutions  for  the  same  airfoil 
and  Mach  number,  but  at  two  different  angles  of  attack,  then  the 
pressure  distribution  for  any  other  angle  of  attack  can  be  found 
by  simple  proportion  and  the  use  of  the  strained  coordinate 
transformation.  Applications  of  this  technique  to  two-dimen¬ 
sional,  multi-parameter  flow  solutions  and  to  three-dimensional, 
multi-shock  problems  are  reported  in  references  1  and  2.  These 
examples  are  concerned  with  the  physical  perturbation  quantities, 
such  as  angle  of  attack,  profile  geometry  changes,  etc. 

More  recently  the  concept  of  "mathematical  perturbations"  has 
been  discussed  (ref.  3)  in  which  the  perturbation  quantities 
characterize  the  difference  between  solutions  of  a  subset  of  the 
full  potential  equation,  namely  the  transonic  small-disturbance 
(or  TSD  equation),  and  the  potential  equation  itself.  The  results 
of  this  type  of  perturbation  solution  indicate  that  such  a  general 
type  of  correction  for  various  subsets  of  equations  is  feasible. 

For  the  present  work  the  perturbation  theory  is  concerned 
with  solutions  of  the  Navier-Stokes  equations.  In  particular, 
three  specific  problems  are  investigated.  In  two  dimensions  the 
application  of  the  physical  perturbation  theory  to  the  Navier- 
Stokes  equations  is  investigated.  Also  the  application  of  the 
mathematical  theory  to  solutions  of  the  transonic  small-distur¬ 
bance  theory  and  the  Navier-Stokes  equations  are  considered.  The 
extension  of  the  mathematical  perturbation  theory  to  three-dimen¬ 
sional  flow  is  also  considered.  In  all  cases  the  strained 
coordinate  technique^  is  employed  to  treat  shock  waves  or  rapid 
gradients  in  the  solutions. 

First  consider  the  physical  perturbations  of  the  Navier- 
Stokes  equations.  The  first  step  in  deriving  a  perturbation 


equation  is  to  write  the  governing  equations  and  their  associated 
boundary  conditions  in  such  a  form  that  only  the  basic  equations 
or  the  boundary  conditions  contain  the  perturbation  parameter, 
otherwise  the  correct  parameter  to  use  may  not  be  obvious.  In 
the  previous  work  for  the  TSD  equation  (ref.  2)  the  perturbation 
parameter  is  contained  only  in  the  boundary  conditions.  For  the 
Navier-Stokes  equations  the  perturbation  parameter  appears  only 
in  the  set  of  equations,  if  these  are  written  in  a  general  body 
conforming  curvilinear  coordinate  system  and  the  variables  non- 
dimensionalized  with  respect  to  free-st  'earn  quantities.  Both 
attached  flows  and  flows  with  shock  induced  separation  are  con¬ 
sidered. 

The  basic  physical  principle  behind  the  "physical"  perturba¬ 
tion  theory  is  that,  relative  to  certain  physical  features  of  the 
problems,  for  example  the  geometric  confines  of  the  airfoil  and 
the  shock  location,  the  physical  processes  throughout  the  pertur¬ 
bation  are  in  some  sense  similar.  For  instance,  if  two  solutions 
are  known  (the  base  and  calibration  solutions)  and  have  shock 
induced  separation  then  an  interpolated  solution  will  also  have 
shock  induced  separation.  Examples  with  attached  flow  and  with 

separated  flow  are  calculated  with  satisfactory  results.  However, 

1  2 

in  contrast  to  the  earlier  potential  equation  work  '  ,  it  is 
found  in  the  present  case  that  the  base  and  calibration  solutions 
should  not  be  too  close,  otherwise  numerical  inaccuracies  in  the 
Navier-Stokes  solution  can  dominate  the  perturbation  quantities. 

The  second  type  of  perturbation  problem  considered  here  is 
the  mathematical  perturbation  correction  technique.  In  the  pres¬ 
ent  case  the  basic  aim  is  to  correct  transonic  small-disturbance 
theory  to  give  results  typical  of  solutions  of  the  Navier-Stokes 
equations.  As  described  fully  in  reference  3,  the  idea  is  to 
determine  a  correction  between  the  two  types  of  equations  by 
considering  the  numerical  results  for  both  equations  for  some 
case.  This  correction  is  then  applied  to  the  TSD  solution  of 
other  near-by  cases.  A  nearby  case  is  assumed  to  mean  an  example 
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with  a  similar  pressure  distribution  to  the  correction  example. 

The  basic  assumptions  behind  this  idea  are,  firstly,  that  a  flow 
is  characterized  by  its  surface  pressure  distribution  and, 
secondly,  that  those  physical  features  modeled  by  the  Navier- 
Stokes  equations  and  not  by  the  TSD  equations  (e.g. ,  viscous 
effects)  for  the  correction  example  are  also  present  in  the 
near-by  cases.  To  a  large  extent  the  validity  of  the  assumption 
is  to  be  determined  by  numerical  experiment. 

A  basic  feature  of  all  perturbation  methods  is  that  the 
perturbation  variables  and  their  derivatives  are  of  the  same 
order  of  magnitude.  In  the  course  of  investigation  of  the  cor¬ 
rection  technique  it  was  found  that  the  shock  movement  due  to  a 
change  in  some  parameter,  e.g.,  angle  of  attack,  was  much  greater 
for  the  TSD  equation  than  for  either  Navier-Stokes  equations  or 
the  Euler  equation.  This  is  due  to  the  larger  shock  strength 
associated  with  the  TSD  equation.  Consequently  an  additional 
perturbation,  within  the  confines  of  small-disturbance  theory, 
is  added  such  that  the  difference  in  shock  strength  between, 
say,  the  Euler  and  TSD  equations  is  the  same  as  for  the  correction 
case. 

The  third  part  of  the  present  work  is  the  extension  of  the 
mathematical  perturbation  correction  technique  (ref.  3)  for 
inviscid  three-dimensional  flow.  To  a  large  extent  this  is  a 
straightforward  development  of  the  work  reported  in  reference  3. 
Both  corrections  for  grid  size,  the  use  of  the  transonic  small- 
disturbance  equation,  and  a  combination  of  both  are  investigated. 

2.  BASIC  PRINCIPLES  OF  THE  PERTURBATION  THEORY 


2.1  General  Concepts 

It  is  usually  assumed  in  perturbation  theory  that  the  form 
of  the  perturbation  parameter  characterizing  the  disturbance  is 
known  or  can  be  easily  found.  For  example,  in  the  earlier  work 


(refs.  1  and  2)  on  the  perturbation  theory,  the  transonic  small- 
disturbance  equation  was  written  in  an  invariant  form  and  the 
"natural"  perturbation  parameters  were  easily  discernible  from 
a  study  of  the  boundary  conditions.  An  example  of  such  a 
"natural"  parameter  is  that,  for  changes  in  Mach  number, 
the  parameter  is 


(1  -  M^) 


where  "A"  denotes  an  increment.  One  problem  that  arose  in  pre¬ 
vious  work  (ref.  4)  was  to  determine  the  precise  nature  of  a  Mach 
number  perturbation  for  the  full  potential  equation,  since  this 
equation  cannot  be  written  in  a  Mach  number  invariant  form.  In 
practice,  the  small-disturbance  perturbation  parameter,  given 
above,  was  used  for  Mach  number  variations. 

In  the  present  work,  the  strained  coordinate  perturbation 
theory  is  to  be  extended  to  treat  the  Navier-Stokes  equations. 

It  is  unlikely  that  these  equations  can  be  written  in  a  form 
independent  of  the  freestream  conditions  and  consequently,  as 
in  the  case  of  the  earlier  full  potential  equation  problems, 
the  choice  of  the  "natural"  perturbation  parameter  for  a  pertur¬ 
bation  in  freestream  quantities  may  not  be  obvious.  Accordingly, 
it  is  proposed  to  write  the  Navier-Stokes  equations  in  a  form 
where  the  boundary  conditions  are  invariant,  and  then  examine  the 
transformed  equations  to  determine  the  correct  choice  of  pertur¬ 
bation  parameters.  A  further  simplification  is  to  use  the  thin- 
layer  (ref.  5)  approximation  to  the  Navier-Stokes  equations. 

This  letter  assumption  is  consistent,  since  these  are  the  equi¬ 
valent  equations  solved  in  most  computer  codes,  whether  explicitly 
coded  or  not. 

The  steady  thin- layer  Navier-Stokes  equations  for  two  dimen¬ 
sions  can  be  written  (ref.  6)  in  a  general  curvilinear  coordinate 
system  as 


M  ^  =  -1  ^ 

9  ?  9  n  e  9  ri 

where  is  the  Reynolds  number,  (^,n)  are  the  computational 
coordinates  which  are  functions  of  Cartesian  coordinates  (x,y) 
with  a  transformation  Jacobian  J,  given  by 


The  transformation  metrics  are  given  by 


C  =  Jy 


h  =  JXr 

y  5 


The  vectors  E,  F,  S  are  defined  as  follows: 


E  =  J 


puU  +  CjjP 

pvU  +  CyP 
(pe  +  p)U 


PMV  +  rijjP 

pvV  +  HyP 
(pe  +  p)V 


'"x  + ''J* “n  1  "x<V„  + Vn’ 

s  =  u(n^  +  nJ)Y^  +  H  ny(n^u^  +  nyV^)  (5 

+  ^  (n  u  +  T)  v)  (n^u  +n  V  ) 

3  X  'y  '  'x  ri  y  n 

where  U,V  are  contravariant  velocities  related  to  the  Cartesian 
velocities  u,v  by 


U  = 

V  =  n.  u  +  n.  V 


(6 


-7- 


p,  p,  e  are  the  density,  pressure  and  internal  energy,  respec¬ 
tively.  is  the  Prandtl  number,  k  is  the  conductivity,  and  a 

is  the  speed  of  sound  of  the  gas.  The  pressure  is  related  to 
e,  p,  u,  V  by 

p  =  (Y-1)  pe  -  +  v^)  (7) 

where  Y  is  the  ratio  of  specific  heats. 

The  computational  coordinates  C/H  ai’e  chosen  such  that  the 
surface  n(x,y)  =  0  represents  the  body  surface.  The  boundary 
conditions  are  the  tangency  and  no-slip  conditions  and  there¬ 
fore. 


V  =  0 

U  =  0  on  T)  =  0 


(8) 


The  far-field  boundary  conditions  are  that  free-stream  conditions, 
denoted  by  the  subscript  prevail. 


If  the  reference  quantities  are  taken  to  be  free-stream 
conditions,  then  the  basic  equations,  eq.  (1) ,  can  be  nondimen- 
sionalized  to  give 


where 


M  +  1£  -  ^ 

an  e  3n 


r 


pU 


puU  + 
pvU  + 
[pe  + 


YM 

oo 

c  p 

y  ym? 


p(Y-l) ]U 


(9) 


(10a) 
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puv  +  n^p 

oo 

—r,  .  -  1 

pvv  +  n„p  — o 

y  YM^ 


(lOb) 


[pi  +  p(Y-l>]V 


y(Ti^  +  n^)u  +  ^n  (nu  +riv) 

X  y  ri  3  'x'  'x  n  y  n 

s  =  j  ^  y(ri^  +  n^)v  +  ^  n  (n  u  +n  v  ) 

X  y  n  3  y  X  n  y  n 

,  2^  2,  v/v  2^  2,  ,-2.-2. 

kp  y (n  +  n  )  ^  Y (Y-i)M  (n  +  n  )  (u  +  v  ) 

r  X  y  8n  2  '  oo '  'x  y  H 

+  4Y(Y-l)M^(n  u  +  n  v)(n  u  +p  v  ) 

I  3  CO  X  y  X  n  y  n 


(lOc) 


where  u  =  u/u„,  v  =  v/u„,  V  =  V/U^,  U  =  U/U„,  p  =  p/p«,  p  =  P/P„f 
e  ==  e/e^,  a  =  a/a^f  and  use  has  been  made  of  the  perfect  gas 
equation 

p/p  =  RT 


and  the  sonic  velocity  relation 


a"  =  YRT 


The  pressure  relation,  eq.  (7) ,  can  be  written  in  the  form 


P  =  (Y-1) 


pe  1 
Y-1  ■  2 


-  -?  -22 
p(u^  +  V  )M^ 


The  boundary  conditions  for  the  set  of  equations,  eqs.  (9)  , 
(10) ,  (11) ,  are  the  tangency  and  no-slip  conditions,  eq.  (8) , 
together  with  the  far-field  conditions 


p  =  p  =  e  =  u  =  v=  l 


It  can  be  seen  then  that,  with  the  exception  of  changes  in 
the  basic  equation  set  and  its  boundary  conditions  are  inde¬ 
pendent  of  the  free-stream  conditions.  Changes  in  geometry  are 

transmitted  through  the  changes  in  the  metrics  ^  /n„/n  and 

X  y  X  y 

the  Jacobian  J. 

The  system  of  equations,  eqs.  (9) -(11),  can  be  written  in 
the  form 

_____  2 

L(p,p,e,u, v,M^,m)  =  0  (13) 


where  m  denotes  the  metric  terms.  Now  consider  a  perturbation 
about  some  state  or  geometry  denoted  by  a  subscript  o,  and  that 
the  perturbation  quantities,  denoted  by  the  subscript  1,  are 
characterized  by  the  parameter  .  First  consider  changes  in 
geometry.  An  expansion  of  eq.  (13)  for  m  then  gives 


o  o' 
3m 


+  eL[pj^,Pj^,eT  ,Ut  ,Vt  ,M^,mJ  =  0 


where  L[  ]  is  a  linear  operator,  related  to  L[  ]  and  which 
depends  on  the  zeroth  order  quantities.  It  can  be  seen  from 
eq.  (14)  that  the  natural  choice  of  perturbation  parameter  is 
Am,  i.e.,  the  change  in  the  geometry  characteristics.  If  the 
general  curvilinear  coordinate  system  is  related  by  a  set  of 
linear  differential  equations  to  the  Cartesian  variables,  then 
Am  is  simply  a  change  in  angle  of  attack,  say,  or  a  change  in 
thickness  ratio.  Hence,  if  two  solutions  are  known  for  two 
cases  that  differ  by  a  one  parameter  variation  then  the  solution 
for  any  other  value  of  this  parameter  can  be  found  by  simple 
proportion. 

A  similar  analysis  to  the  above  but  for  Mach  number  varia¬ 
tions  indicates  that  the  correct  parameter  to  use  for  Mach 
number  changes  is  AM^.  This  is  in  contrast  to  the  earlier 
small-disturbance  result  in  Section  2,  above. 
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Transonic  small-disturbance  theory  can  be  regarded  as  a 
subset  of  the  Navier-Stokes  equations  and  a  question  arises 
regarding  the  relationship  of  these  two  different  Mach  number 
perturbation  parameters,  since  for  certain  flows  the  small 
disturbance  theory  should  be  a  close  approximation  to  the 
Navier-Stokes  equations.  In  the  type  of  interpolation  used 
in  previous  work,  and  also  in  the  present  work,  the  perturba¬ 
tion  parameter  e  appears  in  a  combination,  c/e^,  where  is 
some  known  value  of  e.  In  Table  1  the  effect  of  changing  the 
parameter  for  different  Mach  number  is  shown.  Two  different 
base  and  calibration  Mach  numbers  (which  gives  e^)  are  used 
in  the  computations,  namely  0.8  and  0.83,  and  0.8  and  0.85. 

It  can  be  seen  that,  while  the  values  of  e  for  both  theories 
differ  greatly,  the  ratio  e/e^  for  both  forms  of  e  is  in  good 
agreement.  This  indicates  that  the  small-disturbance  pertur¬ 
bation  is  not  inconsistent  with  the  present  work. 

2.2  The  Strained  Coordinate  Method 

The  interpolation  procedure  outlined  above  is  only  valid 
for  smooth  functions.  A  device  for  treating  discontinuous  func¬ 
tions  using  a  strained  coordinate  system  has  been  derived  in 
previous  papers  (refs.  1  and  3)  and  a  brief  outline  is  given 
below.  Because  of  the  applicability  to  section  4  of  this 
report,  a  more  complete  version  of  the  idea  is  given  here  than 
is  warranted  by  the  simple  interpolation  used  in  this  example. 

The  perturbation  method  outlined  in  the  previous  section 
requires  that  the  changes  due  to  the  correction  are  small.  In 
cases  containing  discontinuities  which  can  alter  location  during 
the  perturbation,  the  changes  are  not  small  in  the  region 
traversed  by  the  discontinuity.  A  means  of  treating  the  prob¬ 
lem  of  perturbations  in  discontinuous  transonic  flow  has  been 
described  by  Nixon  ^ref .  1) ;  an  outline  is  given  below. 


Briefly,  the  idea  is  that  the  problem  is  reformulated  in  a 
strained  coordinate  system  in  which  the  discontinuity  remains  at 
the  same  location  throughout  the  perturbation  and  hence  the  dif¬ 
ficulties  associated  with  moving  discontinuities  do  not  arise 
explicitly.  The  required  straining  is  then  found  as  part  of  the 
solution.  The  basic  equations  in  this  strained  coordinate  system 
are  then  perturbed  about  some  known  solution  to  give  a  linear 
equation  for  the  perturbation  quantities  similar  to  those  dis¬ 
cussed  in  the  previous  section.  Once  the  solution  of  the  linear 
perturbed  equation  is  known,  the  total  perturbed  solution  in  the 
physical  coordinates  is  then  obtained.  The  major  restriction  is 
that  the  discontinuities  must  not  be  lost  or  generated  during 
the  perturbation. 

The  technique  described  above  was  originally  developed  to 
treat  the  discontinuities  which  can  invalidate  a  perturbation 
analysis.  However,  the  technique  can  also  be  applied  to  increase 
the  range  of  application  of  a  valid  perturbation  analysis.  An 
example  from  transonic  airfoil  theory  concerns  the  pressure  dis¬ 
tribution  around  an  airfoil  when  shock  waves  are  present;  such  a 
pressure  distribution  is  sketched  in  figure  1.  The  solid  and 
dashed  lines  denote  two  nearby  solutions  for  the  pressure  dis¬ 
tribution.  The  solution  shock  waves  are  captured,  that  is,  the 
expected  discontinuity  is  smeared  over  a  few  mesh  spacings  and 
denoted  by  CD  and  CD’.  The  method  of  strained  coordinates,  as 
given  in  reference  1,  would  strain  the  x-coordinate  such  that 
the  midpoints  of  the  shock  capture  regions  CD  and  CD'  coincide. 
The  actual  details  of  the  shock  capture  region  are  not  considered 
since  they  are  in  any  event  artificial  phenomena. 

It  can  be  seen  from  figure  1  that  in  the  leading  edge  region 
the  rapid  change  in  the  pressure  distribution  can  cause  large 
pressure  changes  for  a  small  perturbation  if  the  location  of  the 
pressure  rise  shifts  slightly  in  the  x-direction.  This  large 
effect  seriously  limits  the  range  of  validity  of  the  perturbation 
analysis  since  all  pressure  changes  are  assumed  small.  A  method 


of  avoiding  this  difficulty  is  to  strain  the  coordinates  such  that 
representative  point  on  the  AB  and  A'B'  curves  coincide.  This 
then  increases  the  range  of  validity  in  a  similar  way  as  the 
treatment  of  the  shock  waves  in  reference  1. 

A  further  point  concerns  the  treatment  of  the  shock  capture 
regions  CD  and  CD'.  In  the  earlier  applications  (refs.  1  and  2) 
of  the  theory  the  same  mesh  and  differential  equations  were  used 
for  computing  the  pressure  distribution  in  all  examples  and  hence 
the  shock  capture  characteristics  were  essentially  the  same  for 
all  cases.  For  other  problems/  for  example,  those  discussed 
later  in  Section  4,  the  shock  capture  characteristics  may  differ 
substantially  and  it  is  desirable  to  correct  this  behavior. 
Accordingly,  the  coordinates  can  be  strained  such  that  both  the 
points  C,C'  and  D,D'  (the  extremities  of  the  shock  capture) 
coincide.  As  before,  the  actual  flow  details  in  the  region  CD 
are  considered  irrelevant  because  of  the  artificial  nature  of 
the  shock  capture. 

A  more  general  statement  of  the  above  technique  is  as  follows 

1.  If  a  true  discontinuity  is  present,  the  coordinate 
straining  is  such  that  the  location  of  the  discontinuities 
coincide. 

2.  If  there  is  a  shock  capturing  type  of  phenomena,  then 
the  straining  is  such  that  the  extremities  of  the  capture  region 
coincide. 

3.  If  large  gradients  are  present  in  the  solution,  then 
the  coordinate  straining  is  chosen  so  that  a  representative 
point  in  the  region  of  the  large  gradient  coincides. 

These  conditions  constitute  perhaps  a  large  number  of 
requirements  for  the  choice  of  straining.  However,  a  piecewise 
straining  is  perfectly  feasible  provided  the  end  points  of  the 
straining  (which  do  not  move)  lie  in  regions  of  the  solution 
for  which  a  small  perturbation  analysis  is  valid,  for  example, 
in  the  region  BC  of  figure  1. 
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3.  APPLICATION  OF  THE  STRAINED  COORDINATE  SYSTEM 
TO  THE  NAVIER-STOKES  EQUATIONS 


3.1  Theory  for  the  Navier-Stokes  Equations 
Consider  the  set  of  equations,  eq.  (9) . 

3E(g>n)  .  3F(r,,n)  ^  „-i  3s(g,n)  (15) 

3ri  *^e  3n 

Let  both  the  independent  variables  (CiR)  be  strained  such  that 
the  location  of  one  or  more  discontinuities  or  rapid  gradients  are 
held  invariant.  The  new  coordinates  are  denoted  by  (C'fH')  where 


C  =  ^'  +  e6C?3^(C',n’)  +  ... 
n  =  n '  +  eSnrij^  (e '  f  n ' )  +  ... 


(16) 


where  e  is  the  perturbation  parameter,  65, 6n  are  measures  of  the 
movement  of  the  straining  points  and  )  t  are 

(fairly  arbitrary)  straining  functions. 


Now  let  the  dependent  variables  be  expanded  in  a  series  in 
E .  Thus 

E(5,n)  +  cEj^(5',n')  +  ... 


F(5,n)  =  +  eFj^(5’,n’)  + 


(17) 


S(5,n)  =  +  ts^{K\r)')  +  ... 

Using  eqs.  (16)  and  (17)  the  perturbation  equation  for  eq.  (15) 
(the  coefficient  of  e)  is 


3E^ 

W‘ 


3E 

df  1 

3F, 

■3F^ 

3E^ 

65 

O  f. 

35’  H  , 

3n’  . 

n’j 

-  6n 

3n’  %. 

V  * 

35'  ^1  , 

n  J 

==  R 


-1 


3Si  3S^  3S^ 

37^  “  sTT"  '^i  ,  "  FT" 


(18) 
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Now  since  e  is  not  contained  in  eq.  (18)  it  follows  that  a  solution 
to  this  equation  can  be  applied  to  any  value  of  e.  As  in  the  pre¬ 
vious  work  (ref.  2)  the  perturbation  equation  is  not  solved, 
instead  the  solution  of  the  perturbation  equation  is  assumed  to 
be  the  difference  between  two  complete  solutions,  differing  only 
by  a  perturbation  characterized  by  the  parameter  e^.  Hence,  for 
example,  if 


K  =  ,n') 

n  =  n'  +  (?'  ,n' ) 


(19) 


and,  if  for  the  velocity  component,  u(C,ri)/  the  two  solutions  are 
denoted  by  u  (C,ri)  (base  solution)  and  u.  (5,ti)  (calibration  solu- 
tion,  e  =  following  the  previous  work'*, 

u(5,n)  =  u^(5',n*)  +  —  [u^(r,n)  -  UQ(5',n'))  (20) 

and 

e  =  C  +  e  Ci(C’,n*)  =  C  +  ^  (C  -  C') 

°  (21) 

n  =  n'  +  e  6n  Tij^(C’,ri')  =  g'  +  ^  (iT  -  n') 


The  coordinates  ^,ri  are  found  from  eq.  (19)  by  making  the  strain¬ 
ing  points,  e.g.,  C,C'  in  figure  1,  coincide.  A  similar  relation 
to  eq.  (20)  applies  for  all  the  dependent  variables,  the  most 
useful  of  which  is  the  pressure  coefficient  Cp(Cfn).  Thus 

c„(C,n)  =  c  (C’,n’)  +  ^  [c^  (f,n)  -  c  (c',ti’)1  (22) 

P  Po  ^o  Pi  Pi 


In  the  present  work  only  values  on  the  airfoil  surface  (g  =  0)  are 
considered  and  only  the  5  coordinate  is  strained.  The  coordinate 
straining  function  5j^(^*,0)  is  given  by  the  following  equation  or 
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by  a  piecewise  application  of  the  following  equation/  depending 
on  the  number  of  straining  points  required. 


5^{C,0)  =  55^ 


'c'  (l-C)  (C'-C^) 

-4-  * 

+ 

'S' (1-5') (5’-5^)’ 

L^A'i-qxeA-tB'J 

5'  (!-?■')  (E'-E'T 

L.  J 

0  i  X'  j<  1 


=  0 


1  < 


X 


<  0 


This  straining  relation  keeps  invariant  the  leading  and  trailing 
edges  and  two  straining  points,  A,B. 

A  further  result  can  be  obtained  from  eq.  (22)  for  the  lift 
and  pitching  moment  coefficients  ,C  .  Consider  the  lift  coef- 

Lt  lu 

ficient  Cj^,  where 

1 

/  ACp(C,0)dC  (23) 

o 

where  0^(5, 0)  denotes  the  pressure  jump  between  upper  and  lower 
surfaces.  Using  eqs.  (19),  (21),  and  (22)  the  lift  coefficient 
is  then  given  by 


=  1  ( 


(?• ,0')  +  IAC„  (^,0)  -  AC 


vy 

X  l^d?'  +  ^  [df  -  dr/ij 

Now  the  basis  of  the  perturbation  theory  implies  that 


(24) 


[C  (C/n)  -  C  ~  0  (e) 

IdT  -  d^'l  ~  0  (e) 
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and  also  that  d^  =  d^'  +  ^(e) .  Using  these  results  in  eq.  (24) 
it  can  be  shown  that 


C.  = 


AC  (C',0)dr/  +  ^  AC  a,o)dl  - 


—  AC  (^•,0)dC’ 
o  ^o 


(25) 


or 


f-  ^  ' 

o  o  1  o 


(26) 


where  Cl^^  are  the  lift  coefficients  for  the  base  and  calibra¬ 

tion  solutions,  respectively.  A  similar  result  applies  for  the 
pitching  moment  C^^^,  namely 


m 


m 


ic, 


m. 


m 

o 


(27) 


3.2  Discussion  of  Results 

The  variation  of  lift  coefficient  for  a  Korn  airfoil  at 
=  0.755,  for  various  angles  of  attack,  is  shown  in  figure  2. 

Two  separate  cases  are  considered,  namely  flows  with  shock  waves 
and  shock  free  flov;.  For  shocked  flow  the  base  and  calibration 
solutions  were  taken  at  a  =  1.2®  and  a  =  2.66®,  and  for  shock 
free  flow  at  a  =  -0.5®  and  a  =  0.4®.  In  both  cases  the  values 
predicted  by  the  formula  eq.  (26)  agree  satisfactorily  with 
directly  computed  results.  In  figure  3  the  variation  of  C,.  with 
is  shown,  with  base  and  calibration  Mach  numbers  of  0.8  and 
0.85,  respectively.  The  airfoil  is  a  NACA  64A410  section  at  -3.0® 
angle  of  attack.  Again  it  can  be  seen  that  the  simple  formula 
agrees  satisfactorily  with  the  direct  calculation,  except  for 
=  0.75.  The  reason  for  this  is  that  for  all  the  other  Mach 
numbers  the  flow  separates  after  the  shock  wave,  and  this  phenom¬ 
ena  is  implicitly  included  in  the  formula  via  the  and  Cl^.  At 
M^  =  0.75  the  flow  is  attached  and  hence  the  flow  does  not  change 
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smoothly  from  the  base  and  calibration  flows.  A  smooth  change  is 
necessary  for  the  application  of  the  present  theory. 

The  pressure  distribution  around  the  upper  surface  of  two  dif¬ 
ferent  airfoils  calculated  by  the  present  method  are  shown  in 
figures  4-7.  The  first  of  these  airfoils  is  the  NACA64A410  sec¬ 
tion  at  -3.0°  angle  of  attack  and  the  theory  is  applied  for 
different  free-stream  Mach  numbers-  The  base  and  calibration 
Mach  numbers  are  0.8  and  0.88,  respectively.  It  can  be  seen  from 
figures  4  and  5  that  the  present  method  agrees  fairly  well  with 
direct  calculations.  At  the  higher  Mach  number,  M^  =  0.9  the 
agreement  at  the  shock  is  not  quite  satisfactory;  this  may  be 
due  to  the  direct  solutions  changing  character,  either  physically 
or  numerically,  from  the  flow  at  the  base  and  calibration  Mach 
numbers.  These  flows  presented  in  these  figures  have  shock 
induced  separation. 

In  figures  6  and  7  the  pressure  distributions  around  the 
upper  surface  of  the  Korn  airfoil  at  M^  =  0.755  at  different 
angles  of  attack  are  shown.  The  base  and  calibration  angles  of 
attack  are  a  =  1.2°  and  a  =  2.66°,  respectively.  The  agreement 
between  the  present  results  and  the  results  of  direct  calculations 
is  again  satisfactory. 

It  should  be  pointed  out  that  in  the  present  investigation  it 
is  necessary  to  take  the  base  and  calibration  solutions  further 
apart  than  is  the  case  for  potential  flows.  This  is  because  the 
Navier-Stokes  solutions  are  either  not  converged  to  the  same 
level  of  accuracy,  or  that  there  is  some  physical  unsteadiness 
in  the  solution.  Both  these  phenomena  lead  to  "fluctuations'* 
which  can  dominate  the  theoretically  steady  perturbation. 

4.  THE  CORRECTION  THEORY  FOR  NUMERICAL  CALCULATIONS 

The  next  two  sections  deal  with  the  idea  of  correcting 
approximate  solutions  to  give  solutions  typical  of  a  more  accurate 
theory.  It  is  assumed  that  the  more  approximate  solution  can  be 
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generated  much  more  rapidly,  and  less  expensively  than  an  accurate 
solution.  Also,  it  is  assumed  that  the  addition  of  a  correction 
to  the  approximate  theory  does  not  significantly  increase  the 
computing  time.  In  the  next  section  the  idea  of  correcting 
transonic  small-disturbance  solutions  to  give  results  typical  of 
Navier-Stokes  solutions  is  investigated.  In  section  6  the  exten¬ 
sion  of  the  work  of  reference  3,  that  is  the  correction  of  coarse 
grid  solutions  to  fine  grids  and  of  small-disturbance  theory  to 
full  potential  theory  are  examined.  The  basic  idea  behind  the 

3 

correction  theory  has  appeared  elsewhere  but  is  repeated  here 
for  convenience. 


Approximate  equations  are  frequently  in  some  sense  a  logical 
approximation  (rather  than  an  empirical  approximation)  of  an  exact 
equation.  This  means  that  the  approximate  equation  can  be  formally 
obtained  from  the  exact  equation  by  an  expansion  in  some  perturba¬ 
tion  parameter;  powers  of  this  parameter  greater  than  some  specified 
power  are  negligible.  Approximate  equations  can  be  derived  by  such 
an  expansion  are  referred  to  below  as  consistent  approximations. 

The  difference  between  solutions  of  this  approximate  equation  and 
solutions  of  the  exact  equation  is  the  error.  One  example  of  such 
a  consistent  approximation  is  given  by  the  central  finite  difference 
to  (fi  .  Thus 

XX 


4.  -  2(!>  +  4, 

4)  =  - = -  +  Ax  R  , 

2Ax'^ 


(28) 


where  R  is  the  remainder  of  the  truncated  Taylor  series  leading 
to  the  difference  formula.  Finite-difference  representations  of 
a  differential  equation  can  be  represented  by  an  equivalent  dif¬ 
ferential  equation  which  smoothly  approaches  the  exact  differential 
equation  as  the  mesh  size  is  reduced  to  zero.  For  a  finite  mesh 
size  the  difference  between  the  equivalent  and  exact  differential 
equations  can  be  represented  by  a  perturbation  equation  with  a 
perturbation  parameter  characterizing  the  truncation  error. 


In  general,  if  an  approximate  equation  is  a  consistent  approxi¬ 
mation  of  an  exact  equation  then  the  difference  betv/een  the  equa¬ 
tions  can  be  represented  to  a  specified  accuracy  by  a  linear 
perturbation  equation  with  a  perturbation  parameter,  e,  which  is 
a  measure  of  the  error.  If  this  perturbation  equation  is  known, 
it  can  be  solved  using  standard  methods  to  give  a  correction  to 
the  solution  of  the  approximate  equation.  The  main  difficulty  is 
in  actually  determining  the  correct  form  of  the  perturbation 
equation. 

It  is  obvious  that  if  the  exact  equation  can  be  solved  then 
it  is  not  necessary  to  know  the  form  of  the  perturbation  equation 
to  obtain  the  correction  since  the  difference  between  the  exact 
solution  and  the  approximate  solution  gives  the  perturbation  solu¬ 
tion  with  the  addition  of  (formally)  negligible  terms.  If  the 
exact  solution  is  known,  of  course,  then  there  is  little  point  in 
solving  for  the  perturbation  correction.  However,  the  perturbation 
equation  gives  a  formal  measure  of  error  e  between  approximate  and 
exact  equations  and  is  assumed  to  be  small.  If  a  second  perturba¬ 
tion  equation  that  differs  from  the  proper  perturbation  equation 
by  an  error  A  then  the  formal  error  in  using  the  second  equation 
to  compute  the  correction  is  ^  (cA) .  A  problem  which  has  an  error 
described  by  this  second  perturbation  equation  is  referred  to  as  a 
"nearby”  problem.  Its  essential  characteristics  differ  from  the 
real  problem  by  a  small  amount  characterized  by  the  error  A. 

If 

L(<ti)  =  0  is  the  exact  equation 
Lo(i{>o)  =  0  is  the  approximate  equation 
^1^*^!^  =  0  is  the  perturbation  equation 

then 

<1>  =  <|)q  +  +  O(G^)  .  (29) 

If  a  superscript  A  denotes  the  nearby  problem  then,  if 
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L^(<)>^)  =  0 

L^{(j)^)  =  0 
o 

=  0 

Then,  as  before, 

<t>^  =  <^^  +  +  O(e^)  (30) 

Subtraction  of  eqs.  (31)  and  (30)  then  gives,  to  first  order  in  e, 

<1>  =  <|)q  +  (e4>j_  -  e^O^)  +  (9*  -  .  (31) 

If  the  nearby,  approximate,  solution  9^  is  such  that 

U<^1  -  ,  (32a) 

Then  a  good  approximation  to  the  exact  solution  is  given  by 

(()  =  (j)©  +  (<^^  -  4>q)  . 

Since  a  basic  premise  of  perturbation  theory  assumes  that  derivatives 
of  a  quantity  are  the  same  order  of  magnitude  as  the  quantity,  an 
alternative  to  eq.  (32a)  is 

le9i  -  G^<l>n  I  «  |e4),  |  (32b) 

X  X  X 

The  "nearby"  solution  is  defined  by  a  solution  that  satisfies  eqs. 
(30,32).  In  particular  applications  it  is  assumed  that  the  adequacy 
of  a  nearby  solution  can  be  determined  by  examination  of  its 
approximate  solution  in  a  part  of  the  solution  domain,  for  example, 
an  airfoil  surface  rather  than  the  entire  flow  field.  This  assumes 
that  the  characteristics  of  the  complete  solution  can  be  determined 
by  examining  only  a  dominant  part.  A  discussion  on  the  existence 
of  the  nearby  solution  is  given  in  Appendix  B. 
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The  theory  just  outlined  does  not  formally  take  account  of 
shock  waves  or  adequately  treat  solutions  with  rapid  gradients. 

If  either  shock  waves  or  rapid  gradients  are  present  in  the  solu¬ 
tion  then  the  strained  coordinate  method  outlined  in  section  2 
should  be  used.  In  this  case,  each  of  the  solutions  discussed 
above  should  be  written  in  the  appropriate  strained  coordinates 
which  constrain  rapid  gradients  to  the  same  location. 

5.  CORRECTION  PROCEDURES  OF  THE  TRANSONIC 
SMALL-DISTURBANCE  THEORY 

5.1  Theory 

The  transonic  small-disturbance  equation  (TSD)  can  be  written 
in  terms  of  a  perturbation  velocity  potential  4>(x,y)  as 

(1  -  M^)4>  +  d>  =  (y  +  6  (33) 

«■  ^xx  ^yy  ’  ’  “^x’^xx 

where  the  perturbation  velocity  components  in  the  x  and  y  direc¬ 
tions  are  u,v,  respectively  and  are  given  by 

u(x,y)  =  (x,y) ;  v(x,y)  =  (x,y)  (34) 

The  parameter  ”q"  is  a  transonic  parameter  and  within  the  confines 
of  small-disturbance  theory  is  arbitrary.  If  the  slope  of  the 
airfoil  is  given  by  yg(x,±0)  where  "+"  denotes  the  upper  surface 
and  sign  denotes  the  lower  surface  then  the  thin  airfoil 

boundary  conditions  are 

v(x,±0)  =  y^(x,±0)  (35) 

Before  proceeding  with  the  correction  analysis  it  is  perhaps 
pertinent  to  mention  a  few  properties  of  TSD  theory.  The  TSD 
equation  can  be  regarded  as  an  approximation  to  the  conservation 
of  mass  equation;  it  does  not  really  represent  any  of  the  other 
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conservation  laws  such  as  energy,  or  momentum.  The  parameter  q  can 
be  derived  in  a  number  of  ways.  If  the  actual  full  potential  equa¬ 
tion  (ref.  7)  is  truncated  using  the  limiting  idea  of  TSD  theory 

then  q  =  2.  If  the  actual  continuity  equation  is  truncated  (ref.  8) 

2  .  _  2  2 

then  the  term  (Y+1)M  is  replaced  by  (3  +  Y-2M^)M^.  A  commonly-used 
value  for  q  is  1.75,  which  is  Krupp  scaling,  and  is  chosen  to  make 
TSD  solutions  agree  with  full  potential  theory  solutions.  All  of 
these  factors  reduce  to  the  same  value,  (Y+1) ,  in  the  transonic 
limit  of  -i-  1.  The  main  point  of  the  brief  discussion  of  tran¬ 
sonic  small-disturbance  theory  is  that  the  theory  is  flexible  and 
can  be  altered  within  a  certain  range  without  altering  the  formal 
accuracy. 

It  is  interesting  to  briefly  examine  the  behavior  of  TSD  solu¬ 
tions  compared  with  Navier-Stokes  or  Euler  equation  solutions  and 
with  solutions  of  the  full  potential  theory.  The  pressure  distribu¬ 
tions  over  the  upper  surface  of  a  Korn  airfoil  at  =  0.755  and  at 
an  angle  of  attack  of  1.2®  is  shown  in  figure  8.  It  can  be  seen 
that  both  the  TSD  and  full  potential  equation  results  give  a  shock 
location  much  further  aft  than  the  Euler  and  Navier-Stokes  equation 
results. 

The  difference  in  shock  location  betv/een  the  full  potential 
equation  result  and  the  Euler  equation  result  is  about  80%  of  the 
difference  between  the  potential  location  and  the  Navier-Stokes 
location.  Hence,  it  can  be  argued  that  the  cause  of  disagreement 
in  the  shock  location  is  not  primarily  a  viscous  phenomenon,  but 
is  due  to  the  assumption  of  irrotationality  used  in  the  potential 
formulation.  The  reason  for  the  error  in  shock  location  for  the 
TSD  theory  may  be  readily  obtained  from  an  examination  of  the  shock 
strengths  for  the  TSD  equation  and  the  Euler  equation  for  a  given 
pre-shock  pressure.  The  normal  shock  strength  variation  is  plotted 
in  figure  9  and  it  can  be  seen  that,  as  the  magnitude  of  the  pre¬ 
shock  pressure  coefficient  increases,  the  TSD  strength  increasingly 
becomes  much  larger  than  the  corresponding  Euler  strength.  This 
greater  shock  strength  then  requires  the  shock  on  an  airfoil  to 
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move  further  aft  in  order  to  accommodate  it,  in  comparison  to  the 
weaker  Euler  equation  shock.  The  Navier-Stokes  shock  strength  is 
similar  to  the  Euler  strength.  If  a  base  criterion,  that  is,  some 
Mach  number  and  angle  of  attack  are  defined,  and  if  these  criteria 
are  changed  slightly,  it  is  possible  that  the  change  in  shock 
location  for  the  TSD  solution  will  be  considerably  larger  than 
the  change  in  an  Euler  equation  solution.  This  implies  that, 
while  at  the  reference  station  the  difference  in  shock  location 
between  Euler  and  TSD  solutions  may  be  small,  the  difference  in 
the  rate  of  change  of  shock  location  for  both  solutions  is  not 
small.  This  type  of  behavior  can  raise  problems  in  a  perturbation 
analysis,  such  as  the  present  work,  since  perturbation  theory 
usually  requires  that  not  only  are  certain  quantities  small,  but 
that  their  rates  of  change  are  of  comparable  magnitude.  Accordingly, 
this  problem  has  been  addressed  in  the  present  work. 

The  basic  idea  of  the  correction  theory  has  been  outlined 
earlier  in  section  4,  In  the  present  case  the  solution  to  be 
corrected  is  a  TSD  solution  and  correction  is  obtained  from  a 
Navier-Stokes  solution.  From  the  above  discussion  it  is  apparent 
that  if  the  TSD  solution  to  be  corrected  differs  significantly  in 
shock  strength  from  the  "near-by"  solution  used  to  compute  the 
correction,  then  the  new  shock  location  will  be  considerably  in 
error,  even  when  the  correction  for  viscous  effects  is  used.  As 
noted  earlier,  the  TSD  equation  has  a  built-in  flexibility  in 
which  the  transonic  parameter  can  be  varied.  It  is  possible  to 
choose  this  parameter  to  alter  the  shock  strength  or  any  other 
flow  variable  to  agree  with  some  datum  value.  It  is  assumed  here 
that  the  error  in  shock  location  between  a  TSD  solution  and  a 
Navier-Stokes  solution  is  solely  due  to  the  different  shock 
strengths.  Accordingly,  the  TSD  parameter  is  chosen  such  that, 
for  a  given  pre-shock  pressure,  the  difference  betvvreen  the  TSD 
shock  strength  and  the  Euler  equation  shock  strength  is  the  same 
as  in  the  near-by  case  used  to  compute  the  correction.  The  means 
of  accomplishing  this  device  is  as  follows. 
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For  the  near-by  problem  let  the  subscript  o  denote  a  TSD  solu¬ 
tion  for  some  value  of  the  transonic  parameter  q,  and  the  subscript 
1  denote  another  TSD  solution  for  a  second  value  of  the  transonic 
parameter.  The  subscript  2  denote  a  Navier-Stokes  solution  of  the 
near-by  problem.  Finally,  a  subscript  3  denotes  the  TSD  solution 
to  be  corrected.  If  a  denotes  the  shock  strength  this  application 
of  the  usual  form  of  the  correction  theory  will  give  a  new  shock 
strength  where 


where 
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This  may  not  give  the  correct  shock  strength,  and  hence  loca¬ 
tion,  and  so  a  modified  equation  for  the  shock  strength,  within 
the  realm  of  small-disturbance  theory,  is  used,  namely 


0  =  0^  +  ^(^2  “  ^o^ 

The  parameter  e  is  found  by  ensuring  that 

0  -  0^  -  gX(Oj^  -  o^)  =  ^^2  ~  *^o 


Since  -  0^  is  small  and  since  the  "near-by"  solution  is  by 
definition  close  to  the  solution  to  be  corrected  this  device  is 
equivalent  to  tuning  the  small-disturbance  theory  to  get  approxi¬ 
mately  the  correct  shock  jump.  In  eq.  (39)  the  shock  strength  o 
is  determined  at  Cp^,  the  TSD  pre-shock  pressure  coefficient,  by 
the  relation 


o  = 


2Y 

Y+1 


^3 


{M+)^  -  1 


(40) 


where  Mq  is  the  free-stream  Mach  number  and  is  the  pre-shock 

D 

Mach  number  given  by 
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where 


and 


X  =  X2  +  X2-x'  +  e  -  x’ ) 


Xj^  =  x‘  +  6Xg  Xj^ (x* ) 
=  x'  +  6x  X,  (x’ ) 

4.  S2  X 

X,  =  x'  +  6x  X, (x’ ) 
-i  S  ^  ^ 


(44) 


(45) 


where  6X3 6X32/  ^^33  amounts  by  which  the  shock  wave 

moves  for  each  of  the  solutions,  relative  to  the  base  solution 
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(M+)2  = 
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A  similar  relation  holds  for  a^.  The  TSD  strength  is  given  by 


a  =  -2(c;!;  -  C*  ) 

O  P  D 

O 


(42) 


where 


_2  -  <) 


Y+1 


A  similar  relation  holds  for 

Having  obtained  e,  the  final  value  for  the  pressure  coeffi¬ 
cient  on  the  airfoil  surface  is  given  by 


C  (x)  =C  (x,)  +  A[C  (x„)-C  (x’)] 

p  P3  P2  Po 


+  eA[C  (x,)  -  C  (x’)]  (43) 

Pi  1  Po 


where 


and 


X  =  +  X2  -  x'  +  e  (x^  -  x' ) 

X,  =  x'  +  6x  X,  (x' ) 

1  ^ 

x_  =  x'  +  6x^  X, (x* ) 

2  S2  1 

X-  =  x'  +  6x  X, (x’ ) 

3  S3  1 


(44) 


(45) 


where  6X3 5xs2»  <5X33  are  the  amounts  by  which  the  shoc)c  wave 
moves  for  each  of  the  solutions,  relative  to  the  base  solution 
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Cp^{x').  (Note;  The  above  idea  of  a  "self  tuning"  transonic  theory 
is  of  use  in  its  own  right.  A  short  paper  on  this  subject  is 
attached  as  Appendix  C.) 

5.2  Discussion  of  Results 

Several  examples  of  the  correction  technique  have  been  com¬ 
puted  using  both  Euler  equations  and  the  transonic  small  distur¬ 
bance  equation  as  a  base.  With  the  exception  of  the  result  shown 
in  figure  10,  when  the  Navier-Stokes  solver  of  Steger  (ref.  6)  is 
used,  the  Navier-Stokes  and  Euler  equation  solver  used  is  that  of 
Deiwert  (refs.  9  and  10) .  An  earlier  example  for  a  NACA64A410 
airfoil  at  a  =  0  and  =  0.8  is  shown  again  in  figure  10.  The 
basic  calculation  is  an  Euler  equation  solution  and  is  corrected 
for  viscous  effects  from  a  near-by  solution  at  the  sane  angle  of 
attack  and  =  0.82.  A  similar  correction  for  a  Korn  airfoil  at 
=  0.755  and  a  ~  2.66°  is  shown  in  figure  11,  where  the  near-by 
solution  is  at  the  same  Mach  number  and  a  =  1.2°.  Both  these 
examples  have  attached  flow  in  the  near-by  solution  but  the  second 
example  has  mildly  separated  flow  in  the  test  case.  It  can  be 
inferred  from  these  examples  that  the  correction  technique  does 
work  when  Euler  equation  solutions  are  corrected.  It  is  also 
apparent  that  the  error  between  the  Euler  equation  solution  and 
the  Navier-Stokes  equation  solution  is  not  great  and  that  these 
examples  are  perhaps  not  a  good  test  of  the  method.  Furthermore, 
the  computation  time  of  the  Euler  solver  is  about  half  that  of  the 
Navier-Stokes  solver  and  the  prospective  overall  savings  in  compu¬ 
tation  cost  is  not  as  dramatic  as  would  be  the  case  of  correcting 
TSD  equation  solutions.  Accordingly,  most  of  the  effort  was  put 
into  investigating  the  correction  procedure  for  the  TSD  equations. 

In  figure  12  the  pressure  distribution  around  a  Korn  airfoil 
at  M  =  0.755  and  a  =  1.2°  is  shown.  The  base  result  is  a  TSD 

CO 

calculation  using  Spreiter  scaling  and  the  correction  is  computed 
for  the  same  airfoil,  at  the  same  Mach  number,  with  a  =  1.7°.  In 
figure  13  the  pressure  distribution  around  the  Korn  airfoil  at 
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=  0.755  at  a  =  2.66°  is  shown.  The  near-by  solution  is  the  same 
as  in  the  previous  example.  It  may  be  seen  that  for  both  these 
examples  that  the  corrected  result  is  fairly  satisfactory  although 
the  shock  location  is  not  quite  right.  However,  the  results  are 
a  considerable  improvement  on  the  uncorrected  TSD  solution  using 
Spreiter  scaling  and  a  greater  improvement  over  those  computed  with 
Krupp  scaling.  It  should  be  noted  here  that  the  present  method 
effectively  chooses  its  own  transonic  scaling  through  the  use  of 
eq.  (35) .  It  is  interesting  to  investigate  what  effect  the  use 
of  eq.  (35)  has  on  the  solution.  Hence,  the  preceding  case  was 
computed  without  the  use  of  the  shock  strength  scaling  correction 
and  the  result  shown  in  figure  14.  It  can  be  seen  that  the  shock 
strength  correction  is  indeed  an  essential  part  of  the  theory. 

The  previous  example  is  one  primarily  of  attached  flow.  An 
attempt  is  made  here  to  try  to  compute  a  separated  flow  example. 
Although  separation  is  really  an  unsteady  process,  the  results  used 
here  are  computed  using  a  steady  version  of  Deiwert's  Navier-Stokes 
solver.  The  flow  around  an  18%  thick  biconvex  airfoil  at  a  =  0° 
and  =  0.775  is  shown  in  figure  15.  The  near-by  solution  is  the 
same  airfoil  at  the  same  angle  of  attack  but  at  =  0.79.  The 
agreement  between  the  corrected  result  and  the  direct  result  is 
fairly  good  although  the  pressure  in  the  separation  region  is  too 
high.  However,  the  corrected  result  is  again  a  considerably 
improvement  over  the  TSD  solution.  It  is  worthwhile  noting  that 
the  turbulence  model  had  to  be  the  same  both  for  the  near-by  case 
and  for  the  direct  comparison  computation;  different  turbulence 
models  can  produce  considerably  different  results. 

It  is  one  of  the  implications  of  the  present  method  that  the 
near-by  solution  need  not  be  for  the  same  airfoil  on  the  test  case 
only  that  its  basic  characteristics,  i.e.,  rapid  gradients  should 
be  similar.  Accordingly,  an  example  for  a  NACA0012  airfoil  at 
a  =  0°  and  =  0.8  was  computed  using  the  correction  obtained  for 
the  upper  surface  of  the  Korn  airfoil  at  =  0.755  and  a  =  1.2°. 
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The  resulting  pressure  distribution,  shown  in  figure  16  agrees 

surprisingly  well  with  the  direct  result.  The  NACA0012  result 

is  mildly  unsteady  and  the  range  of  unsteadiness  is  indicated  by 

the  symbols  on  the  figure.  An  attempt  was  made  to  use  the 

correction  obtained  for  a  NACA64A410  airfoil  at  M  =0.8  and 

00 

a  =  -3.0°  for  the  biconvex  airfoil  discussed  above.  This  result 
is  shown  in  figure  17  and  it  can  be  seen  that  although  the  use 
of  the  correction  procedure  upgrades  the  TSD  solution,  the  overall 
improvement  is  not  really  satisfactory.  The  main  reason  for  the 
error  in  the  pressure  in  the  separation  region  seems  to  be  that 
the  inviscid  TSD  solution  for  the  NACA64A410  airfoil  does  not  have 
as  large  a  pressure  gradient  at  the  trailing  edge  region  as  the 
biconvex  airfoil  solution.  Incidently,  it  should  be  noted  that 
the  examples  shown  in  figures  12  and  13  are  formally  beyond  the 
bounds  of  small  disturbance  theory. 

In  general  then  it  appears  that  the  correction  theory  works 
for  correcting  Euler  equation  solutions  and  for  TSD  equation  solu¬ 
tions  for  attached  flow.  The  shock  locations  can  be  in  error  and 
this  may  be  due  to  the  extreme  nature  of  some  of  the  examples.  It 
also  appears  that  the  use  of  a  different  airfoil  for  the  correction 
computation  can  work  although  in  extreme  cases  like  that  shown  in 
figure  17  the  improvement  may  be  minimal.  Also,  from  the  one 
example  computed,  it  seems  as  if  massively  separated  flows  can 
be  computed  although  more  work  is  required  to  determine  the  range 
of  applicable  "near-by"  solutions. 

Most  of  the  examples  discussed  above  are  fairly  extreme  in 
their  degree  of  flow  complexity.  One  of  the  difficulties,  in 
abortive  attempts  to  compute  the  flow  over  the  NACA64A410  airfoil 
at  a  =  -3.0°  and  the  Mach  numbers  used  in  figures  4  and  5,  is  that 
the  TSD  solution  gives  a  shock  location  at  the  trailing  edge  even 
for  quite  high  values  of  the  transonic  parameter,  thus  destroying 
any  chance  of  using  the  shock  location  correction  of  eq.  (35) . 

The  overall  conclusion  reached  from  the  small  range  of  examples 
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computed  is  that  the  correction  theory  does  appear  to  work,  but 
that  a  much  greater  range  of  examples  both  less  extreme  than  the 
present  ones  and  with  massive  separation  need  to  be  investigated. 

A  major  difficulty  in  implementing  the  suggestion  economically  is 
to  predict  in  advance  what  kind  of  flow  a  Navier-Stokes  solution 
will  exhibit.  A  further  point  that  should  be  investigated  is 
some  means  (probably  semi-empirical)  of  identifying  vvhen  a 
separated  or  attached  flow  correction  should  be  used. 

6.  THE  CORRECTION  THEORY  IN  THREE-DIMENSIONS 

6.1  General  Comments 

The  extension  of  the  mathematical  perturbation  theory  to  three 
dimensions  is  a  straightforward  development  of  the  theory  described 
in  section  4.  The  correction  for  a  given  problem  is  obtained  by 
computing  one  exact  and  one  approximate  solution  and  taking  the 
difference.  The  strained  coordinate  method  is  used  to  treat  shock 
waves  and  any  rapid  gradients  in  the  flow.  The  only  fundamental 
difference  between  two  and  three  dimensional  applications  of  the 
theory  is  that  in  three  dimensions  provision  has  to  be  made  for  a 
spanwise  straining  to  treat  possible  shock  intersection  points. 

The  examples  considered  in  the  present  investigation  all  con¬ 
cerned  the  ONERA  M6  wing.  It  is  obviously  desirable  to  treat  the 
theory  over  a  wide  range  of  wings,  but  the  relative  unavailability 
of  wing  input  data  and  shortage  of  time  inhibited  the  use  of  wings 
other  than  the  M6  wing  as  test  cases.  However,  the  pressure  dis¬ 
tribution  around  this  wing  at  typical  transonic  speeds  is  fairly 
complex  and  usually  contains  two  shock  waves  intersecting  inboard 
of  the  wing  tip,  together  with  a  large  pressure  gradient  in  the 
vicinity  of  the  leading  edge.  In  the  test  cases  various  combina¬ 
tions  of  angle  of  attack  and  wing  twist  were  used  to  try  to  produce 
a  relatively  wide  range  of  perturbations. 


I 
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An  important  decision  in  the  use  of  the  theory  is  to  determine 
the  number  and  location  of  the  straining  points.  A  typical  section 
pressure  distribution  for  the  M6  wing  is  shown  in  figure  18.  It 
can  be  seen  that  the  leading  edge  shock  wave  is  weak  and  is  con¬ 
siderably  smeared  fay  the  finite-difference  schemes;  this  can  cause 
difficulties  in  identification.  Consequently  it  was  decided  to 
simply  treat  the  leading  edge  region  as  a  region  of  high  gradients 
and  to  pick  some  identifier  (not  necessarily  the  shock)  to  hold 
invariant.  The  ultimate  choice  is  the  location  of  the  bottom  of 
the  first  compression  gradient  (point  A  in  figure  18) .  The  other 
points  held  invariant  in  the  straining  functions  are  the  beginning 
and  end  of  the  capture  region  of  the  strong  aft  shock,  points  B,  C, 
respectively  in  figure  18.  The  leading  and  trailing  edges  are  also 
held  invariant.  In  addition  to  this  section  data  a  shock  inter¬ 
section  point  on  the  span  is  also  held  invariant.  It  should  be 
noted  that  the  solution  accuracy  is  unchanged  if  a  previously 
specified  rapid  gradient  does  not  exist  at  some  straining  point, 
since  the  straining  does  not  degrade  a  normally  valid  or  accurate 
perturbation.  One  problem  that  arises  in  a  diverse  straining, 
such  as  the  one  necessary  in  the  present  problem,  is  that  the  usual 
cubic  straining  function  used  in  section  3  can  give  strained  coordi¬ 
nates  that  tie  off  the  wing.  In  order  to  remove  this  tendency  a 
piecewise  version  of  the  linear  straining  g'ven  in  reference  2  is 
used.  Fo.-  a  chord  length  of  unity,  the  three  straining  points  A, 

B,  C,  the  straining  function  is 

C^(C,0)  =  C'  +  C'<SC;^H(C^-C’)  +  [6?;g(r-C^)  - 

X  H(^'-r,^)H(^-C')  +  I65^(C*-C^)  - 

X  H(?'-e^)H(f.^-C')  +  6C^(l-^')/(l-F,^)H(C'-C<l:)  (46) 
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where  H{^)  is  the  step  function,  and  6^',  6^',  and  6^'  are  the 

A  15  C 

perturbations  in  the  points  5',  5'. 

A  15  \»i> 

This  straining  function  is  used  in  the  computation  of  all  the 
test  cases. 


6.2  Discussion  of  Results 

The  correction  for  all  the  examples  is  taken  from  various 
computations  of  the  ONERA  M6  wing  at  =  0.84  and  a  =  3.06“ 
and  zero  twist  angle.  Both  full  potential  calculations,  using 
Jameson's  finite  volume  method  (ref.  11),  and  modified  small- 
disturbance  calculations,  using  a  variant  of  the  Bailey-Ballhaus 
(ref.  12)  code  are  performed.  Calculations  are  carried  out 
relatively  coarse  and  fine  grids.  The  type  and  magnitudes  of 
errors  involved  in  the  various  approximations  are  indicated  in 
figure  19.  The  span  stations  n  are  such  that  n  =  0  denotes  the 
wing  root  and  0  =  1  denotes  the  %ving  tip.  As  pointed  out  in 
section  4,  the  main  criteria  for  both  correction  results  and  the 
approximate  or  coarse  grid  results  is  that  they  must  contain  all 
of  the  essential  features  of  the  flow.  For  example,  in  problems 
of  the  type  sketched  in  figure  18,  a  coarse  grid  solution  must 
capture  the  leading  edge  compression  in  an  adequate  fashion  other¬ 
wise  the  location  of  the  forward  straining  point  cannot  be  found. 
Another  interesting  variation  of  the  criteria  for  validity  is  that 
the  convergence  level  for  all  solutions  on  the  same  grid  must  be 
the  same.  If  the  convergence  level  (and  hence  the  truncation  error) 
is  different  for  different  cases,  then  the  equivalent  perturbation 
equation  for  the  correction  is  not  the  same,  and  the  present 
correction  technique  is  invalid.  It  is  possible  however  to 
envisage  a  correction  which  would  treat  different  convergence 
levels. 

The  first  two  examples,  shown  in  figures  20  and  21,  indicate 
the  application  of  the  correction  technique  for  grid  size.  The 
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approximate  solutions  are  obtained  using  the  finite-volume  method 
(ref.  11)  and  are  on  a  grid  of  dimensions  80x8x16;  the  fine  grid 
has  dimensions  of  160x16x32.  The  results  for  =  0.84  and  a  =  3.06 
are  shown  in  figure  20  for  various  span  stations  n  and  it  can  be 
seen  that  the  agreement  of  the  present  method  with  direct  calcula¬ 
tions  is  fairly  good.  In  figure  21  the  pressure  distribution 
around  the  upper  surface  of  the  M6  wing  at  =  0.84,  a  =  3.06 
and  with  3.0°  of  wing  twist.  Again  it  can  be  seen  that  the  agree¬ 
ment  between  the  corrected  solutions  and  the  directly  calculated 
solutions  is  adequate. 

In  figure  22  and  figure  23  the  corrections  for  the  use  of  the 
TSD  equation  is  shown.  The  TSD  solution  uses  the  grid  embedding 
technique  of  Boppe.  In  figure  22  the  pressure  distribution  for 
=  0.84  and  a  =  3.06  is  shown,  and  in  figure  23  the  pressure 
distribution  for  =  0.84,  a  ==  4.0°  and  1.5°  of  twist  is  shown. 

In  both  cases  the  agreement  with  direct  calculations  is  adequate. 

Finally,  in  figure  24  and  figure  25,  the  correction  for  grid 
size  and  the  use  of  the  TSD  equation  is  investigated.  The  TSD 
solutions  are  computed  on  a  grid  with  31  streamwise  points  and  the 
full  potential  equation  on  the  grid  of  160  16  32  as  before.  In 
figure  24  the  pressure  distribution  for  =  0.84,  a  =  3.06  is 
shown  and  it  can  be  seen  that  the  agreement  with  the  direct  calcu¬ 
lation  is  not  really  adequate  away  from  the  wing  root.  This  is 
almost  certainly  due  to  the  inability  of  the  coarse  grid  TSD 
solution  to  adequately  represent  a  main  feature  of  the  final 
solution,  in  this  case  the  "bulge"  in  the  pressure  distribution 
in  the  region  of  the  leading  edge  recompression.  It  should  be 
noted,  however,  that  the  present  technique  still  upgrades  the 
solution,  especially  the  shock  location.  In  figure  25  the  pressure 
distribution  for  a  =  4.0°  and  1.5°  of  twist  is  shown.  It  can  be 
seen  that  these  corrected  results  compare  fairly  well  with  the 
results  of  direct  calculations. 

In  summary  then,  the  correction  technique  derived  for  two- 
dimensional  flows  in  earlier  reports  does  appear  to  work  in 


j 

A 


three-dimensions  although  care  must  be  exercised  in  the  choice  of 
the  grid  size.  A  major  criterion  of  the  work  is  that  only  one 
wing,  the  ONERA  M6  wing,  is  used  in  the  tests.  The  reason  is 
given  elsewhere,  but  the  method  should  be  more  comprehensively 
tested  for  other  wings.  A  second  point  is  that  the  computer  code 
used  in  the  present  research  is  not  entirely  automatic.  A  future 
code  should  have  an  automatic  identification  method  for  the 
straining  points. 


7 .  CONCLUDING  REMARKS 

The  main  object  of  the  present  work  is  to  extend  some  recent 
development  in  perturbation  theories  of  transonic  flow.  One  topic 
is  the  extension  of  the  physical  perturbation  theory  to  treat 
solutions  of  the  Navier-Stokes  equation.  A  second  topic  is  the 
investigation  into  the  possibility  of  correcting  lower  grade 
inviscid  solutions  for  viscous  effects  using  the  mathematical 
perturbation  theory.  Thirdly,  the  extension  of  the  mathematical 
perturbation  theory  for  potential  flows  to  three  dimensions  is 
considered.  The  extension  of  the  physical  perturbation  theory  to 
solutions  of  the  Navier-Stokes  equations  is  straightforward,  the 
only  additional  fact  to  appear  is  that  the  necessary  base  and 
calibration  solutions  should  not  be  too  close  together,  otherwise 
the  perturbation  quantities  can  be  seriously  degraded  by  the 
numerical  error  in  the  solution. 

Application  of  the  correction  technique  to  upgrade  potential 
equation  solutions  to  include  viscous  effects  introduced  some 
interesting  problems.  The  main  difficulty  proved  to  be  due  to 
the  fact  that,  for  relatively  strong  shock  waves,  the  potential 
theories  give  shock  locations  that  are  much  further  aft  than  those 
predicted  by  the  Euler  or  Navier-Stokes  equations.  This  raised 
problems  in  the  application  of  the  concept  of  near-by  solutions, 
since  the  sensitivity  of  some  of  the  examples  to  small  geometry 
changes  could  allow  a  flow  change  from  one  with  a  relatively  weak 
shock  to  one  with  a  strong  shock.  This  led  to  the  necessity  for 
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a  correction  for  the  purely  inviscid  errors  due  to  the  use  of  the 
potential  equation  rather  than  the  Euler  equations.  Examples 
computed  with  the  subsequent  correction  theory  compare  adequately 
with  direct  calculations. 

The  extension  of  the  mathematical  correction  theory  to  three- 
dimensions  is  a  straightforward  development  of  the  t\TO-dimensional 
theory.  Within  certain  limits  the  perturbation  method  appears 
satisfactory.  However,  in  all  of  the  topics  discussed  above,  it 
is  desirable  that  a  more  comprehensive  testing  procedure  be  per¬ 
formed  over  a  wider  range  of  airfoils  and  wings  in  order  to  assess 
the  range  of  applicability  of  the  theories. 
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TABLE  1.-  VARIATION  OF  e  AND  —  WITH  MACH  NUMBER 
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Figure  2.-  Variation  of  with  a,  Korn  airfoil 

=  0.755. 


Figure  3.-  Variation  of  Cl  with  M 
NACA64A410  airfoil,  a  =  -3.0°. 
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— “Direct  calculation 
O  Perturbation  method 


Figure  4.-  Pressure  distribution  around  the  upper 
surface  of  a  NACA64A410;  =  0.85,  a  =  -3.0. 
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Figure  7.-  Pressure  distribution 
surface  of  a  Korn  airfoil; 


around  the  upper 
=  0.75,  a  =  2°. 
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figure  10,-  Pressure  distribution  around  MACA64A410  airfoil 
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Figure 


11^ _  Pressure  distribution  around  the  upper  surface 
of  a  Korn  airfoil;  =  0.755,  a  =  2.66, 
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Navier-Stokes  solution 
Uncorrected  TSD  solution 


Figure  12.-  Pressure  distribution  around  the  upper  surface 


of  a  Korn  airfoil:  M  =  0 

^  CO 


755,  a  =  1. 2' 


Navier-Stokes  solution 


Figure  15.-  Pressure  distribution.- around  a  18-percent 
Biconvex  airfoil:  M  =  0.775,  a  =  0°. 
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Figure  20.-  Continued. 
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Figure  21.-  Pressure  distribution  over  the  upper  surface 
of  the  ONERA  176  wing;  M  =  0.84,  a  =  3.06,  ^twist  ~ 
(grid  correction,  fu°Cl  potential  equation)  . 
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Figure  21.-  Continued. 
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21.-  Continued. 


Figure  22.-  Pressure  distribution  around  the  upper  surface 
of  ONERA  M6  wing;  =  0.84,  a  =  3.06° 

(TSD/full  potential  equation  correction) . 
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Figure  23.-  Pressure  distribution  over  the  upper  surface 
of  the  ONERA  M6  wing;  a  =  4°,  ^twist  “  1.5° 
(TSD/full  potential  equation  correction) 
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Figure  23.-  Continued. 
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23.-  Concluded. 
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Figure  24.-  Concluded. 
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Figure  25.-  Continued. 
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Figure  25.-  Concluded. 
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APPENDIX  A 

OBSERVATIONS  ON  THE  STPJilNED  COORDINATE 
METHOD  FOR  TRANSONIC  FLOWS 


* 

By  David  Nixon 

Nielsen  Engineering  &  Research,  Inc. 
Mountain  View,  California 


Introduction 

1  2 

Recently  papers  '  concerning  a  strained  coordinate  method  for 

discontinuous  transonic  flow  problems  have  appeared  in  the 

literature.  In  particular,  the  application  of  this  technique  for  a 

linear  interpolation  or  extrapolation  between  two  known  results 

has  raised  some  questions  regarding  the  relationship  of  the  strained 

coordinate  technique  to  normal  interpolation  procedures.  Some 

confusion  also  seems  to  have  occurred  regarding  the  treatment  of 

3 

the  Osv/atitsch -Zierep  "shock  foot  singularity."  It  is  the  purpose 
of  this  note  to  examine  these  questions  and  (hopefully)  finally 
resolve  them. 


Analysis 

In  figure  1  a  typical  velocity  distribution,  for  transonic 
flow,  u(x),  is  sketched.  This  is  the  curve  ABCD.  A  shock 
wave  is  present,  denoted  by  the  line  BC,  with  the  shock  foot  singu¬ 
larity  at  the  point  C.  At  C  the  derivative  ^  is  logarithmically 
singular.  A  second  transonic  flov/  related  by  some  small  parameter 
change  to  the  first,  is  represented  by  the  curve  A3CD.  '  The 
purpose  of  a  linear  interpolation  procedure  is  to  obtain  a  third 

Index  Categories:  Transonic  Flow;  Analytical  and  N'oraerical  Methods. 
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solution,  A  B"  C"  D  from  the  first  tv;o  (known)  solutions. 

Let  the  difference  between  the  AfiCD  and  AB '  C '  D  solution  be 

characterized  by  the  small  parameter  e  and  the  difference  bet'ween 

o 

A3"C"D  and  ABCD  characterized  by  e.  The  shock  locations  C,  C',  C" 
are  denoted  by  x^,  x^,  x^.  In  both  the  strained  coordinate  method 
and  linear  interpolation  the  new  shock  location  is 


X 


s 


(X-X-) 


(1) 


Consider  now  the  velocity  distribution  u(x)  in  the  region  AB  where 

no  discontinuities  or  singularities  appear.  The  strained  coordinate 
2 

method  gives 

u(x)  =  u^(x*)  +  ~  [Uj^(x)  -  u^(x')]  (2) 


where  u^(x')  refers  to  ABCD,  Uj^(x)  refers  to  AB'C'D.  The  physical 
coordinate  is  x'  and 


X  =  x'  +  6x  x,  (x’ ) 

S  X 

(3) 

X  =  x'  +  —  6x  x, (x' ) 
e  si 

(4) 

o 


where 


and 


X  (x‘)  =  -  X') 

^^l'^  ^  X'  (1  -  x;) 


5x  =  X  -  X ' 
s  s  s 


A  normal  linear  interpolation  procedure  gives 

u(x)  =  +  —  [u^{x)  -  Uq(x)] 

o 


(5) 


(6) 


Now  from  equations  (3)  and  (4) 

X  =  x'  +  —  5x  X,  (x' ) 
e  si'' 
o 

X  -  X  -'r  (~  -  1)  6x  Xj^(x') 

o 


and  thus,  by  a  Taylor  series  expansion, 


u^(x)  =  U^(x-)  +  6XgX^(x')  ^  +  ... 

o 


u,  (x)  =  u,  (x)  +  - 1)  6x  X,  (x')  ~  + 

^  ^  ®  3x 


(7) 


Substitution  of  equation  (7)  into  equation  (6)  and  retaining  only 
linear  terms  in  Sk^  and  e  gives 

u(x)  =  u^(x')  +  -^  - 1)  (5XgXj^(x) 

+  ^  [Uj^Cx)  -  u^(xM]  (8) 


i!i  ^ 

9x  " 


If,  as  will  be  the  case 


in  region  AB, 


0  (e) 


(9) 


then  the  second  term  on  the  right  of  equation  (8)  can  be  neglected 
and  it  may  be  seen  that  the  strained  coordinate  method  and  linear 
interpolation  are  equivalent.  A  similar  equivalence  occurs  in  the 


region  C'D  provided  the  point  in  question  is  not  too  close  to  the 

3ui  3u 

singularity  at  C  .  This  v/ould  make - large  compared  v;ith  ^ 

3x 

and  would  invalidate  the  assumption  of  equation  (9) . 

Consider  now  interpolation  or  extrapolation  in  the  region 
BB'.  An  extrapolation  for  u^(x')  can  be  made  using  a  Taylor  series 
expansion  about  some  point  x^  close  to  B.  Thus 


Uo(x) 


Uo(x^)  +  (X 


3u 

x' )  - ^ 

[ax- 


A  similar  analysis  to  that  previously  given  leads  to  the  result 


u{x)  =  u^(x^)  +  - 


O  [“1 


(x)  -  u  (x') 
o  p 


+  f  F--1  vr  -  top  F--1  )  (11 

o  >>  o  9x  ^  o  X 

^  P 

3u^  3u 

If  — =■  and  3—°  ai’e  sufficiently  close  and  if  Ax  ~  0  [6x  x,  (x’)] 

then  the  third  term  in  equation  (11)  is  negligible  and  again  the 
linear  interpolation/extrapolation  is  equivalent  to  the  strained 
coordinate  method.  A  similar  analysis  can  be  applied  to  the  region 


C"C'  in  which  case  if  x'  is  close  to  C 

P 


3u^ 
_  3x- 


is  very  large  and 


interpolation/extrapolation  is  not  equivalent? to  the  strained 

r9u,i 


irdinate  method  since 


would  not  be  sufficiently  close  to  the 


3u  -  -  p 

rr-or^'  regular  3 — in  C"C’.  Indeed,  it  is  difficult  to  see  how 

oX 

r  C"C'  can  be  obtained  by  linear  interpolation/extrapolation 


,011  Id  involve  an  expansion  through  the  singularity  at  C. 


It  is  probably  possible  hoi^ever  to  obtain  a  solution  in  CC  if  the 
shock  jump  relations  are  used  at  B"  to  get  values  at  C”  and  a  hig2ier 
order  interpolation  used  in  C'C". 

Consider  now  the  behavior  of  the  strained  coordinate  method 
at  the  shock  foot.  The  behavior  of  the  velocity  at  C,  C  is  given 
by  the  form 


u^  (x)  =  ^  -  x^)!!,n{x  -  Xg)e 

u  (x')  =  u  (x')  +  a  (x’ -  x*)An(x’-  x') 
o  o  s  o  s  s 


(12) 


where  ^  values  at  B'  and  B  respectively  and 

*^1'  *^o  constants  depending  on  the  shock  strength  and  the  air¬ 

foil  curvature.  Using  the  strained  coordinate  method,  equation  (2), 
gives  for  values  just  downstream  of  C" 

u(x)  =  u  (x'  )  +  a  (x'-  -  x*) 

o  s__  os  s 

^  “l^^  -  -  x^)  (13) 

-  a  (x*-  x')Jln(x'-  x') 
o  s  s 

From  equations  (3)  and  (5)  it  can  be  seen  that 

(x  -  X  )  =  x' -  X '  +  6x  (x,  (x‘)  -  1) 

S  S  S  X 

If  (x'  -  Xg)  =  5  then,  using  equation  (5), 

(X  -  X  )  =  611  +  a(x’)l  (14) 

O 


where 
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Substitution  of  equation  (14)  into  equation  (13)  gives 


u(x)  =  u  (x  )  +  a^6£n5  +  —  (x  )  -  u  (x^  ) 
o  s_  o  [Is  o  i 


+  aj^6[l  +  a(x’)l£n6[l  +  a(x')] 


-  ci  6ZnS  > 

°  J 


(16) 


Since  ox  is  assumed  small 
s 


£n6[l  +  a(x')]  -  £n6 


Then 


u(x)  =  |u^(x 


;)  +  r-  ^ 

o  ~ 


} 


+  Ict^  +  ^  (Oj^ll  +  a(x')]  -  aQ}6^n6 
o 


(17) 


The  form  of  the  singular  part  of  u(x)  at  C"  can  be  found  by  the 
above  analysis  to  be 


+  ^  a(x' )]  fn5  [l 


where  02  is  some  constant. 

Since  a(x')  is  small  this  can  be  reduced  to 

a2[6£n6  +  {6£n6)  ~  a(x')j 


Thus 


6£n6  = 


6[l  + 


e 


a(x’) 


+ 


e 


a(x’) 


Substitution  in  equation  (13)  gives 


u  (x)  = 


+  fa„  +  ^  <“l  -  »o>  +  ^  -  ^]a(x’)  (aj^ 

L  o  o  '■  o-' 

X  (5^1  +  a(x')j£n^|^l  +  ~  a(x*)j^ 


Since  we  are  dealing  in  sraall  perturbations  it  is  obvious  that 

0  0 

(a,  -  a  )  is  small  and  hence  the  terra  —  (1  -  — )a(x')  [o.-.  -  a  )  can 
1  o  e  e  1  o 

o  o 

be  neglected  in  equation  (18) .  It  can  then  be  seen  that  the  strained 
coordinate  system  not  only  transfers  the  shoc)c  to  its  correct  loca¬ 
tion,  from  equation  (1) ,  but  also  introduces  the  correct  shock  foot 

singularity  v/ith  a  strength  [a^  +  ■—  (a^  -  a  )],  obtained  by  linear 

o 

interpolation  of  the  strengths  at  C,  C .  Thus  in  the  strained 
coordinate  method  the  treatment  of  the  shock  foot  singularity  is 
consistent  with  the  treatment  of  the  smoother  parts  of  the  solution. 

Conclusions 

The  strained  coordinate  interpolation  method  used  for  transonic 
flows  has  been  compared  with  normal  interpolation/extrapolation 
procedures.  It  is  found  that  both  methods  are  essentially  equiva¬ 
lent  in  smooth  regions  of  the  solution.  However,  normal  linear 
extrapolation  does  not  appear  to  be  applicable  in  the  region  just 
behind  a  shock  wave  where  the  acceleration  is  infinite.  The 
strained  coordinate  method  does  move  the  shock,  and  its  associated 
shock  foot  singularity  to  the  correct  location  and  scales  the 
strength  of  the  singularity  according  to  linear  interpolation. 
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APPENDIX  B 


COMMENTS  ON  THE  EXISTENCE  OF  A  "NEAR-BY"  SOLUTION 


In  the  preceding  work  the  term  "near-by  solution"  is  used. 
The  definition  is  very  loose,  being  simply  "a  solution  that  has 
all  the  essential  features  of  the  basic  solution."  In  this 
appendix  an  attempt  is  made  to  define  a  near-by  solution  on  a 
quantitative  basis. 

Consider  a  differential  operator,  L(  )  where 

L((!))  =  0  (1) 

Let  an  approximate  representation  of  the  system  described  by 
equation  (1)  be  given  by 

L($)  +  eL3^(^)  =  0  (2) 

where  e  is  a  measure  of  the  difference  between  approximate  and 
exact  solutions.  $  is  the  solution  to  the  approximate  equation 

L2(i)  =  0  (3) 

and  Lj^  (  )  is  a  differential  operator.  If 

([)  =  $  +  £())'  (4) 

where  e  is  a  small  parameter,  then 

L(4))  =  L($)  +  eL((J)')  +  ...  =  0  (5) 

where  L(  )  is  a  linear  differential  operator  whose  coefficients 
may  be  functions  of  $,  and  ({>'  is  the  correction  to  ({>. 
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If  equation  (5)  is  subtracted  from  equation  (2) ,  then  an 
expression  for  the  correction  '  is 

L((f)')  =  L^($)  (6) 

Now  let  there  be  two  solutions,  a  base  solution,  denoted  by  the 
subscript  1,  and  a  near-by  solution,  denoted  by  the  subscript  2. 
Thus,  from  equation  (6) 


L2((|)^)  = 

Subtraction  of  these  equations  give 


(7) 


Now  generally  in  transonic  flow  problems  L(  )  is  a  mixture  of 
first  and  second  derivatives  (no  cross  derivatives) .  Let  L(  ) 
be  given  by 

2  2 

L(  )  E  a($)  +  e($)  +  Y(^)  -h  +  (5) 

The  boundary  condition  for  the  perturbation  equation,  equation 
(6) ,  is  given  by 


({)y  =  B  on  some  boundary  C  (10) 

The  precise  location  of  the  boundary  C  may  not  be  too  critical; 
i.e.,  a  mean  surface  boundary  condition  could  be  used. 
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Equations  (8)  and  (9)  can  be  combined  to  give 


^“2  "  “l^'*’2xx  ®l*‘^2yy  *^2  ”  "'^1^ 


+  («2-«i)fiy  +  q<ii>  -  qi»2> 


with  boundary  conditions 


'♦i-<’2>y  “  ®2  -  ®l 


on  some  mean  of  Cj^,C2 


~  ^2'  i“l 


Ul  -  82!  «  IBjlI 
-  Y2I  <<  Iy^I 
|6^  -  62I  «  I6j 
iB^  -  B2I  «  |B^t 


|Li(<f)i)  -  L^((})2 


)  1  I  «  1  iL^(^i)  1  I 


where  |  [  1  1  is  some  suitable  norm  for  ((^)  ,  then  equations 

(11)  and  (12)  reduce  to 


0,4''  +  3,4''  +  Yi4>'  + 

1  XX  1  xy  'lx  1  y 


4>’  ~  0 

y 


where 


4>'  =  <!>1  -  ^'2 


The  solution  to  equations  (14)  and  (15)  is 


(J)’  =  0 


or 

<t)j  -()>•  =  0 

Formally  this  implies 

1‘t'i  "  ‘^’2'  I'*’!' 

which  in  turn  implies  that  the  correction  used  for 

the  problem  denoted  by  subscript  2  provided  in  equation  (13) 
is  satisfied. 


Example ; 


Consider  the  problem  of  truncation  error  in  solutions  of 
the  TSD  equation 

+  A  =  4  ij)  (1 

XX  yy  ^x  XX 

In  this  case 


„  =  -  4,^ 

B  =  1 


Y  =  <|), 


6  =  0 


L^($)  =  Ax^(B^ 


<i>  ) 


ill 

9x^ 


+  Ay' 


(1 


Hence,  the  correction  ijrocedure  is  valid  if 
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'•4  -  ♦!=:  -  ®2  ^  *2x1 


<< 


1^1  -  ^ 


lx' 


^Ix  ■  1^ 


2x' 


lx' 


-  ‘*>lx> 


4~  4~ 

9  <j),  3  (()- 


9^4>1  2 

-  (B^  - 
9x^  ^ 


9  ^2 


2x' 


9x 


<< 


2  ~ 

-  ‘^ix^  7^ 
9x 


9y 


3y 


<< 

9y^ 

(18) 


The  above  formula  are  dependent  on  large  gradients  either 
vanishing  or  being  constrained  to  the  same  location.  The 
strained  coordinate  method  takes  account  of  the  latter  problem. 


In  practical  applications  it  is  assumed  that  the  error 
bounds,  equation  (13) ,  can  be  deduced  from  values  on  the  air¬ 
foil  surface.  This  does  not  seem  unreasonable. 


i 

1 


1 


APPENDIX  C 


TRANSONIC  SMALL  DISTURBANCE  THEORY  WITH 
STRONG  SHOCK  WAVES 

David  Nixon* 

Nielsen  Engineering  &  Research,  Inc. 

Introduction 

The  most  conunon  methods  of  predicting  aerodynamic  character¬ 
istics  at  transonic  speeds  are  either  the  Transonic  Small 

1  7 

Disturbance  (TSD)  theory  or  the  Full  Potential  Equation  theory 

(FPE) .  The  more  accurate  Euler  equation  solutions  are  expensive 
to  obtain  although  for  flov^7s  with  strong  shock  waves  such  solu¬ 
tions  are  essential.  The  FPE  theory  requires  that  the  flow  is 
irrotational  and  treats  the  wing  boundary  conditions  exactly 
(numerically) .  The  TSD  theory  is  an  approximation  to  the  FPE 
theory.  One  advantage  of  the  TSD  theory  is  the  flexibility  in 
deriving  the  approximate  equation.  This  flexibility  is  generally 
utilized  by  a  choice  of  a  transonic  scale  parameter.  The  basic 
assumption  of  irrotationality  in  both  these  theories  is  only 
valid  when  the  flow  is  shock  free  or  contains  only  weak  shocks. 
Both  TSD  and  FPE  solutions  are  in  satisfactory  agreement  v/ith 
realistic  Euler  equation  solutions  provided  that  the  basic 
restriction  to  weak  shock  waves  is  not  violated.  The  thin 
wing  boundary  conditions  can  also  introduce  errors  into  the  TSD 
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* 

Research  Scientist,  Associate  Fellow  AIAA. 


-2- 


solutions.  If  the  flow  has  strong  shock  waves,  hov/ever,  then 
there  is  considerable  disagreement  among  all  three  theories. 
Generally  the  predicted  shock  locations  for  the  potential  theories 
are  much  further  aft  than  for  the  Euler  equations.  The  problem 
addressed  in  this  paper  is  to  examine  the  error  in  the  shock 
location  in  the  TSD  theory  in  two-dimensional  flow  and  to 
derive  a  correction  procedure  within  the  confines  of  small 
disturbance  theory.  The  basic  hypothesis  of  the  present  theory 
is  that  the  error  in  shock  location  is  primarily  due  to  the 
stronger  shock  strength  predicted  by  TSD  theory  compared  to  that 
of  the  Euler  equations.  The  technique  uses  two  TSD  solu¬ 
tions  with  different  scaling  parameters  and  an  interpolation 
scheme  derived  for  discountinuous  transonic  flows  to  give  a 
corrected  shock  strength. 


Analysis 


The  TSD  equation  for  the  perturbation  velocity  potential, 
4>(X/y)/  at  a  free-stream  Mach  number  M^,  is 


M^)4>  +  *  =  (y  +  1)M*^(1)  6 

oo'  ^xx  ^yy  '  »^x^xx 


where  y  is  the  ratio  of  specific  heats  and  q  is  the  transonic 

4 

scaling  parameter.  The  two  most  commonly  used  values  of  q  are 
2  (Spreiter  scaling)  and  1.75  (Krupp  scaling).  The  pressure 
coefficient  Cp(x,y)  is 


Cp(x,y)  =  -2  ()>j^(x,y) 


(2) 
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Associated  with  Eq.  (1)  are  the  usual  tangency  and  far  field 
boundary  conditions.  The  weak  shock  jump  conditions  for  Eq.  (1) 
are 


a  1  -  IC  -  (y  + 


6  B  +  tan  (J)  [I  D 


(3) 


where  H  D  denotes  a  jump  through  a  shock  wave  and  ^  is  the  angle 
betv/een  the  shock  normal  and  the  x-axis.  For  a  shock  normal  to 
the  free  stream  the  shock  strength,  is  defined  as 


a 


c 

p 


(4) 


where  C^,  are  the  pressure  coefficients  just  ahead  of  and 
behind  the  shock  and 


*  “2(1  -M^) 

C  =  - 

^  (y  + 


(5) 


Consider  now  the  Euler  equation  normal  shock  jump,  in 

+  5 

terms  of  and  which  is  given  by 


"f,  =  <»e  -  » 


lYM, 


2  p 


(6) 
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where  the  upstream  shock  Mach  number  is  given  by 


Eqs.  (4)  and  (6)  are  shown  for  =  0.755  in  figure  1  and 
it  can  be  seen  that  as  increases  the  discrepancy  between 

and  o_  increases.  Note  that  different  transonic  scalings 

X  Jbi 

'fc 

not  only  give  a  different  value  of  C  but  generally  a  different 

P 

value  of  Cp.  Thus,  for  different  scalings  the  shock  strength 
may  vary  considerably. 


The  error  in  the  shock  location  in  the  TSD  solutions  seems 


to  be  primarily  due  to  the  error  in  the  shock  strength  as 
exhibited  in  figure  1.  If  the  TSD  equation  is  altered,  still 
within  it's  formal  accuracy  bounds,  such  that  the  shock  jump 
approximates  the  Euler  equation  shock  jump  then  the  resulting 
equation  is  a  better  compromise  in  representing  the  actual 
flow.  The  reason  for  this  statement  is  that  by  matching  the 
shock  jiomp  the  new  equation  implicitly  introduces  an  additional 
mechanism,  formally  negligible,  that  cancels  the  rotationality 
errors  in  a  potential  formulation. 


If  the  correct  shock  strength  is  known  and  if  two  TSD 
solutions  with  different  scalings  are  also  known  then  a  TSD 
solution  with  correct  shock  strength  may  be  estimated  using  a 


-5- 


linear  combination  of  the  known  solutions.  Thus,  if  and 

o,_  are  the  shock  strengths  of  the  TSD  solutions,  then  a  param- 
eter  e  can  be  found  such  that 


where  Og  is  the  estimated  Euler  equation  shock  strength.  If  e 
is  known  then  the  pressure  distribution  can  be  found  using  the 
strained  coordinate  method^.  Thus, 

C_,(x,y)  =  C  (x’,y)  +  e  IC  (x,y)  -  C  (x',y)]  (9) 

P2 


where  C  (x,y) ,  C  (x,y)  are  the  pressure  coefficients  from  the 

Pi  ^2 

TSD  solutions  and 


X  =  x’  +  6x  X, (x' ) 


X  =  X'  +  e  (x  “  X' ) 


(10) 


\?here  6x  is  the  change  in  shock  location  between  first  and 
o 

second  TSD  solutions.  The  function  Xj^(x’)  is  usually  taken 
to  be 


Xj^(x*) 


where  x^  is  the  shock  location  for  the  first  TSD  solution, 
s 
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ja 


The  main  problem  now  is  to  estimate  the  shock  strength  a^. 
This  is  given  by  Eq.  (6)  if  is  knovi^n.  In  TSD  solutions 
Spreiter  scaling  gives  good  agreement^  for  moderately  strong 
shocks  and  Krupp  scaling  gives  good  agreement  for  weak  shocks. 
Hence  it  seems  reasonable  to  take  to  be  the  average  of  the 
Spreiter  and  Krupp  values  of  the  preshock  pressure  coefficients. 
However,  the  justification  for  this  choice  really  lies  in  the 
accuracy  of  the  final  results.  Hence,  given  o^.  from  Eqs.  (6) 
and  (7)  the  surface  pressures  can  be  found  from  Eqs.  (8)  ,  (9)  , 
and  (10) . 


Results 

The  present  method  was  first  tested  for  a  weak  shock 
example,  namely  the  flow  over  the  NACA  64A006  airfoil  at  zero 
angle  of  attack  and  =  0.875,  to  see  if  the  present  method 
would  give  the  accurate  Krupp  scaling.  This  result  is  compared 
to  an  Euler  equation  result^  in  figure  2  and  it  can  be  seen  that, 
apart  from  the  shock  capture  characteristics,  the  present 
method  does  agree  satisfactorily  with  both  the  Euler  solution 
and  the  Krupp  solutions.  An  example  for  a  flow  with  a  strong 
shock  is  shown  in  figure  3.  In  figure  3  the  flow  around  a 
Korn  airfoil  at  =  0.755  and  1.7“  angle  of  attack  is  compared 

Q 

to  an  Euler  equation  solution  and  it  can  be  seen  that  the 
agreement  of  the  shock  location  predicted  by  the  present  method 
is  satisfactory,  although  there  is  an  error  in  the  leading  edge. 


which  is  almost  certainly  due  to  the  use  of  thin  airfoil  boundary 


n 


conditions  in  the  TSD  theory.  Incidently,  if  Spreiter  scaling 
is  used,  the  TSD  result  is  almost  coincidental  with  the  present 
result. 


Concluding  Remarks 

A  method  to  effectively  choose  a  transonic  scaling  to 
place  shock  waves  computed  by  small  disturbance  theory  at  the 
location  predicted  by  Eulor  equation  solutions  has  been 
developed.  The  technique  does  effectively  correct  the  shock 
location  although  discrepancies  in  the  leading  edge  region 
still  persist.  This  is  probably  due  to  the  use  of  thin  wing 
boundary  conditions. 
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Figure  1.  Variations  of  normal  shock  strength  with  pre-shock 
pressure. 

Figure  2.  Pressure  distribution  around  a  NACA  64A006  airfoil 
a  =  0°,  =  0.875. 

Figure  3.  Pressure  distribution  around  the  upper  surface  of 
Korn  airfoil;  =  0.755,  a  =  1.7®. 
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